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o : 

^ Based on density-functional perturbation theory we have computed the photoelastic tensor of a 

' model of sodium silicate glass of composition (Na20)o. 25(8102)0. 75 (NS3). The model (containig 84 

atoms) is obtained by quenching from the melt in combined classical and Car-Parrinello molecular 
dynamics simulations. The calculated photoelastic coefficients are in good agreement with experi- 
mental data. In particular, the calculation reproduces quantitatively the decrease of the photoelastic 
^\ ' response induced by the insertion of Na, as measured experimentally. The extension to NS3 of a 

phenomenological model developed in a previous work for pure a-Si02 indicates that the modulation 
upon strain of other structural parameters besides the SiOSi angles must be invoked to explain the 
O ' change in the photoelstic response induced by Na. 
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I. INTRODUCTION 



Photoelasticity is of interest from the fundamental point of view as well as for several technological apphcations in op- 
tics and microelectronics. For instance, photoelasticity in a-Si02 is known to cause a reduction of fiber Bragg gratings 
efficiency^ and to produce a loss of resolution of pure silica lenses used in photolytography^. Moreover, photoelastic- 
[ ity is directly related to the Rayleigh scattering coefficient in single-component glasses through the Landau-Placzek 
\ relation, and is thus one of the main contributions to loss in state-of-the-art silica fibers for telecommunications'^. A 
Q ■ systematic experimental study of photoelasticity in pure and modified silica glass has been reported by Schroeder 
O in the early 80s.^ Brillouin scattering measurements have shown that the modification of silica glass with alkali or 
alkali-earth ions (Li, Na, K, Ca, and Mg) reduces sizably the photoelastic coefficients.'' This result is in contrast with 
I ■ the prediction of simple models (Lorenz-Lorentz) based on the observed increase in density and refractive index with 
^ , the modifier concentration. As expected, a related reduction of the Rayleigh scattering coefficient was later measured 
C<1 ' for these glasses,'^ and Lines thus suggested their use as ultra-low loss glasses for telecommunication applications.^ 
\ In a recent work^ we have shown that the photoelastic coefficients of crystalline and amorphous pure Si02 can be 
computed with good accuracy within density functional perturbation theory. A phenomenological model based on ab- 
initio data allowed us to identify the microscopic parameters which rule photoelasticity in pure a-Si02. In the present 
paper we have applied the same ab-initio framework to compute the photoelastic tensor of a sodium silicate glass 
with composition (Na20)o.25(Si02)o.75 (NS3) aiming at identifying which modification cither structural or electronic 
: induced by the insertion of sodium is mostly responsible for the change in photoelastic response. Models of NS3 
d glass containing 84 atoms have been generated by quenching from the melt in combined classical and Car-Parrinello 
^ \ molecular dynamics simulations. The calculated photoelastic coefficients are in good agreement with experimental 
I ' data. In particular, the calculation reproduces quantitatively the changes of the photoelastic response induced by 
■ O the insertion of Na, as measured experimentally. In order to identify the microscopic mechanims responsible for the 
^ : change of the photoelastic response induced by Na, we extended to NS3 the phenomenological model developed for 
pure a-Si02 in ref.^. The paper is organized as follows. In section II we describe our computational framework. In 
section III we report the details of the structural and elastic properties of our model of NS3. The results on the 
calculated photoelastic tensor and their interpretation in terms of a phenomelological model are presented in sections 
IV and V, respectively. Section VI is devoted to discussion and conclusions. 



II. COMPUTATIONAL DETAILS 



A model of sodium silicate glass of composition (Na20)o.25(Si02)o.75 (NS3) has been generated within a combined 
classical and ab-initio framework which has been previously used to generate models of pure a-Si02^'^. Models of the 
glass have been generated by quenching from the melt in classical molecular dynamics simulations and then annealed 
for few ps within ab-initio Car-Parrinello molecular dynamics*. The same method has been used by Ispas et al.^ to 
generate a theoretical model of sodium tetrasilicate glass ((Na20)o.2(Si02)o.8: NS4) achieving good agreement with 
experiments both in the structural and electronic properties. Ispas et al.^ used a modified BKS-potential^*^ extended 



f 



to sodium silicates by Horbach ct al}^. Conversely, we have adopted the empirical potential developed by Ovicdo ct 
al}^ which is an extension to sodium silicate glass of the forceficld introduced by Vashista et al}^ for pure amorphous 
silica. It consists of a short-range two-body interaction, long-range coulomb interactions, and a three-body term which 
enforces the directionality of the Si-0 covalent bond. The Na-Si, Na-0 and Na-Na interactions have been modeled 
by Oviedo et al}"^ by coulomb interactions plus a repulsive short-range term. The effective charges of the coulombic 
potential are: 1.0 for Na, 1.6 for Si, and -0.971 for oxygen which enforces charge neutrality in NS3 (note that there is 
a misprint in the oxygen charges in table I of ref.^^). A model of liquid NS3 has been prepared by inserting 7 Na20 
and 21 Si02 units in a cubic box (a=10.531 A) at the experimental density of glassy NS3 (2.427 g/cve?). The system 
is equilibrated at 6800 K, cooled to 3800 K in 50 ps and then equilibrated at the final temperature for 5 ns. The 
sample is then quenched at room temperature in 5 ns (quenching rate of 7-10^^ K/s) and equilibrated at 300 K for 50 
ps. This model has then been annealed at 600 K in Car-Parrinello simulations for 1.3 ps and then quenched at 300 
K in 0.15 ps. Structural properties have been averaged over a NVE nm at room temperature, 1.1 ps long. The ab- 
initio simulations are based on density functional theory in the local density approximation (LDA)^** as implemented 
in the code CPMD*. Norm-conserving pscudopotential for Si and Na have been used. Non linear core corrections 
are included in Na pseudopotential^^. An ultrasoft pscudopotential^^ has been used for oxygen. Kohn-Sham (KS) 
orbitals are expanded in plane-waves up to a kinetic cutoff of 27 Ry. Integration of the BZ has been restricted to the 
r point. To study the dielectric properties, the structures generated by Car-Parrinello simulations have been then 
optimized with norm-conserving pseudopotentials and a larger cutoff of 70 Ry. We have computed the dielectric and 
photoelastic tensors within density functional perturbation theory DFPT^^, as implemented in the code PWSCF and 
PHONONS^^. The photoelastic tensor is defined by 

Ae,^-^ = Pijkir]ki (1) 

where Sij are the components of the optical dielectric tensor, -qki is the strain tensor. Only the electronic contribution 
is included in Sij, also indicated as e°°. Experimentally, it corresponds to the dielectric response measured for 
frequencies of the applied field much higher than lattice vibrational frequencies, but lower than the frequencies of the 
electronic transitions. The photoelastic coefficients have been calculated by finite differences from the dielectric tensor 
of systems with strains from -1 % to -|-1 %. The components of the photoelastic tensor will be expressed hereafter in 
the compressed Voigt notation. 

The exchange-correlation functionals available in literature (LDA and generalized gradient approximation (GGA)) 
usually underestimate the electronic band gap and overestimate the electronic dielectric tensor up to 10-15 %^®. 
This discrepancy can be corrected scmi-empirically by applying a self-energy correction, also referred to as a scissor 
correction, which consists of a rigid shift of the conduction bands with respect to the valence bands^°. This procedure 
has been used successfully to reproduce the photoelasticity of Si^^, GaAs^^ and quartz^^. However, as shown in 
our previous work^, it turns out that even within simple LDA, the error in the photoelastic coefficients for several 
polymorphs of silica is smaller than what expected on the basis of the error in the dielectric constant itself. The 
scissor correction has thus been neglected in the calculation on NS3 reported here. 



III. STRUCTURAL PROPERTIES 



As pointed out in previous works'^^'^'^^, the insertion of sodium modifies the topology of the network by the 
formation of non-bridging oxygens (NBO). The number of NBO coincides roughly (and in our small cell, exactly) 
with the number of Na atoms. All silicon atoms in our model are 4- fold coordinated and bonded at most with one 
NBO. The presence of NBO drastically reduces the number of medium sized rings (5-8 Si atoms per ring) which are 
usually predominant in the ring-size distribution of pure silica (Fig. 1). 
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FIG. 1. The ring size distribution of (a) a model of pure a-Si02® and of (b) the NS3 model after CPMD annealing. 

The structural properties of the final ab-initio model of NS3 are compared to those of the classical model (before 
ab-initio annealing) in Figs. 2, 3 and 4. Pair correlation functions (g(r)) and radial coordination numbers are shown 
in Fig. 2. The average Si-0 bond for NBO is slightly shorter than for bridging oxygen ions (BO) as shown in Fig. 5 
which report gsio{r) resolved for BO and NBO. 

The ab-initio annealing produces a slight broadening of the pair correlation and angle distribution functions (ADF, 
Fig. 4). In fact, the three-body term in the classical potential assigns a tetrahedral geometry too stiff with respect 
to the ab-initio results. The O-Si-O ADF is still centered aroimd the tetrahedral angle (109.5°), while the Si-O-Si 
ADF is centered at 139.3°, which is smaller than the average Si-O-Si angle in pure silica. Furthermore, in contrast 
with the classical models of pure silica generated with the BKS potential (cfr. ref.^'^'^) no shift in the maximum of 
the g{r) (SiO, 00 and SiSi) is observed upon ab-initio annealing. In fact, the Vashista potential is fitted directly on 
the experimental g{r) of amorphous silica and the agreement with ab-initio bondlengths is better than for the BKS 
potential. However, the ab-initio annealing shifts outwards the peak of the Na-0 correlation function (from 2.4 to 
2.5 A). By separating the contribution of NBO and BO to the pair correlation functions (Fig. 3) we can see that the 
Na atoms are closer to NBO than to BO. On average, a NBO is coordinated with 2.9 Na ions and a BO with 1.3 Na 
ions. On the other hand, a Na ion is on average coordinated with about 2.9 NBO and 3.35 BO (cfr. Fig. 2). 
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FIG. 2. Pair correlation functions of the NS3 model. Bold solid and dashed lines correspond to ab-initio and classical MD 
simulations, respectively. The thin solid line is the ab-initio radial coordination number. 
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FIG. 3. Partial pair correlation functions for NBO and BO in the classical (dashed line) and ab-initio (bold line) models of 
NS3. 



Looking at the Na-Na pair correlation function it is worth noting that the broad peak at 3.5 A present in the classical 
MD model, disappears in the CPMD simulation. Thus, sodium does not tend to segregate or to form clusters, but it 
is homogeneously distributed in the network. 
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FIG. 4. Angular distribution functions. Bold solid and dashed lines correspond to ab-initio and classical MD simulations, 
respectively. 
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FIG. 5. Silicon-oxygen pair correlation function for NBO (dashed line) and BO (solid line) of the ab-initio model of NS3. 

Our results on the structural properties of NS3 are close to those obtained by Lspas et al® for a model of NS4 within 
a similar theoretical framework. 

After the ab-initio annealing performed with the softer Vanderbilt pseudopotentials, the final structure has been 
further optimized with norm-conserving pseudopotentials (for the calculations of the dielectric properties). The cell 

geometry has been optimized at fixed volume allowing orthorhombic distortions of the initially cubic supcrcell such 
as to produce a diagonal stress tensor. The residual anisotropy in the stress (a) for the optimized ratios of cell edges 
b/a= 1.04 and c/a = 1 is 



-502.8 -10.9 
C7= I -10.9 -495.1 
4.03 -18.6 




kbar 



(2) 



The large negative stress in eq. 2 is due to the so-called Pulay stress. The b/a and c/a ratios obtained in this way 
at the experimental equilibrium density have been held fixed and the volume varied to generate the equation of state 
reported in Fig. 6. 
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FIG. 6. Ab-initio equation of state of the NS3 model. 



The calculated E{V) points have been corrected for the discontinuities due to the incomplete basis set following 
the prescription given in rcf.^'' and then fitted by a Murnaghan function^^. The resulting equilibrium density (peq), 
bulk modulus {B) and derivative of the bulk modulus with respect to pressure {B') are peg = 2.468 g/cm'^ (exp. 
2.427 g/cm-^), B = 44.1 GPa, B' = 1.6. The insertion of Na produces a marginal softening of B with respect to pure 
silica^. In fact in our previous work^ on a model of a-Si02 of similar size (81 atoms) we obtained in our previous 
work^ we have obtained Peg = 2.31 g/cm^, B = 46.2 GPa, B' = —6.28. The value of B' , negative in pure silica, turns 
to positive in NS3, which means that the structural response to compression is different in silica and in NS3. Na-0 
interaction is probably responsible for the change in B' . 

The response to strain (7722,^/33) of some average structural parameters of NS3 and of a-Si02 (81-atoms supercell) 
are compared in table I. The change in the SiOSi angles upon strain is substantially smaller in NS3 than in pure 
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a-Si02, nevertheless the response of the Na-0 distances to strain is rather relevant, if compared to the average Si-O 

distance changes, both in pure sihca and in NS3. 



TABLE I. The derivatives of the SiO and NaO bond-lengths, SiOSi angles and Si-Si nearest neighbors vector distances with 
respect to strains r/22 and 7733 in the models of NS3 and ar-Si02®. Pi(Si-Si) denotes the projection of the average Si-Si vector 
distance on the i-th axis. NBO and BO indicate non-bridging and bridging oxygen atoms, respectively. 
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IV. DIELECTRIC AND PHOTOELASTIC PROPERTIES 



The dielectric tensor of the NS3 model optimized at its equilibrium density is: 

/ 2.449 0.008 -0.006 \ 

e= 0.008 2.438 0.031 , (3) 
\ -0.006 0.031 2.474 / 

The average theoretical dielectric constant of NS3 and of pure silica (cfr. Ref.[6]) are compared to experimental data 
in table . 

TABLE IL Theoretical and experimental dielectric constant of NS3 and pure a-Si02. " Ref.[6] . 

DFPT exp.^ 

NS3 2.454 2.2.36 

a-SiO^ 2.292" 2.125 



The increase of s due to sodium insertion is quantitatively reproduced by our calculations. The computation of the 
photoelastic tensor has been performed by finite differences of the dielectric tensor by applying strains of ±1%. 

Results on the photoelastic cocfScients are compared with experiments^ in table III. The pu and pi2 coefficients 
have been obtained by averaging over different components which should be equal in a homogeneous model, i.e. 
pii = (p3;5 +P22)/2 and pi2 = {pi2 +P13 + P23 +P,32)/4. The agreement with experimental data is of the same quality 
as for pure a-Si02. In particular, the calculation reproduces quantitatively the change in photoelastic coefficients 
observed experimentally upon insertion of Na, namely a large decrease of the off-diagonal pi2 and a smaller increase 
of pn. 



TABLE III. Ab-initio photoelastic coefficients of the NS3 model compared with experimental data, with photoelastic coeffi- 

cioiits (if pure silica coiriputod l)y DFPT and with measured ones. 
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V. PHENOMENOLOGICAL MODEL OF PHOTOELASTICITY 



A. Pure silica glass 

In a previous work^, we have developed a phenomenological model of the dielectric properties of silica based on 
ab-initio data. The main features of the model are briefly outlined here. Its extension to NS3 is described in the next 
section. In rcf.® wc have assumed that the dielectric response of silica polymorphs could be embodied in an ionic 
polarizability tensor of the oxygen ions, whose value is assumed to depend on the SiOSi angle only. The dielectric 
susceptibility x (e = 1 + ^ttx) can be obtained from a site dependent oxygen polarizability, aj, as 



(4) 



where V is the volume of the unit cell containing N oxygen ions and B is a 3Nx3N matrix consisting of 3x3 blocks 
Bik defined as 



(5) 



and 
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'''ik - 



ViV, — = 



(6) 



where R are Bravais lattice vectors defined by the shape of the supercell for models of glasses or by the unit cell for 
crystalline phases, rik is the distance between sites i and k in cells separated by R (see Ref.^ for details). 




FIG. 7. Sketch of the Si-O-Si unit and of the bond contributions to polarizability. 

The polarizability of the oxygen ions is a function of the SiOSi angle 9 and can be expressed in terms of the 
polarizability of the SiO bonds as (cfr. ref.^) 

c(6l) +7cos2(6»/2) 

c{e)+-^sin^{9/2) | (7) 

aT>{0) 



with c = 2{a.L + cut + a^O/S, 7 = — ut and a^, and q;t' arc the longitudinal and the two transversal 
polarizabilities of the Si-0 bond in Fig. 7, respectively. The contribution of each polarizable SiOSi unit to the 
dielectric susceptibility is 



R[a{e)R, 



(8) 



where Ri is the rotation matrix that operates the transformation from the local reference system represented in figure 
7 to the absolute reference system of the solid, in which the i-th Si-O-Si unit is embedded. The parameter 7 has been 
assumed independent on 6 and set equal to the value obtained from the fitting of the ab-initio Raman spectrum of 
a-quartz in ref.^^ (7 = 9.86 a.u.). As discussed in our previous work^, the functions c{9) and aT'{0) have been fitted 
on the dielectric properties of a-cristobalite at different densities. The results reported are shown in Fig. 8. 
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FIG. 8. The functions c{9) (solid line and circles) and Qt'(^) (dashed line and squares), which assign the oxygen polarizability 
(see text). The data are the result of the fitting on the dielectric tensor of a-cristobalite at different densities. 

This model allows the calculation of the photoelastic tensor given by the derivative of the dielectric susceptibility 
tensor x (eq- 4) with respect to the strain tensor rjxn as 




rp da 



-^=aRi + Ria--=- 



A + 



(9) 



where the matrix A is defined as -A = {L — B_)~^ (see eq. 4), and the arguments in square brackets indicate 3x3N 
matrices. The change of the oxygen polarizability with strain can be expressed in terms of dc{6)/d9 and dax' (0) / 09 
deduced from Fig. 8. All the other derivatives can be obtained by finite differences. The phenomenological model 
outlined above has been shown to reproduce satisfactorly the dielectric and photoelastic tensors of several silica 
polyphorms (quartz, a-cristobalite, /3-cristobalite, a-Si02) at normal conditions and at high density^. 



B. Extension to sodium silicates 



We can now make use of this phenomenological model for pure a-Si02 to identify which term in Eq. 9 is mostly 
affected by the presence of sodium. As shown in Ref.^ the change in e and in the density alone within a simple 
Lorenz-Lorentz model (cfr. section II) can not account for the change in the photoelastic coefficients measured 
experimentally upon Na insertion. This is inferred by comparing the quantity pde/dp measured experimentally 
{{pds/dp)obs = e^(Pii + 2pi2)/3) with the result of the Lorenz-Lorentz model 

de _ {e-l)ie + 2) 
^^^'"^ 3 

For pure a-SiOs, {pde/dp)obs=l-003 and {pds/dp)LL=l-5i7, while for NS3 {pds/dp)obs=0-937 and 
{pde/dp)LL=^-74:6. Thus, the Lorenz-Lorentz model predicts an increase in the photoelastic response with Na content 
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which is in contrast with the experimental observation. An increase in the dependence of the molecular polarizability 
upon strain {da/drj in Eq. 9), neglected in the Lorentz-Lorenz model, is required to account for a decrease in pi2. 
By still keeping valid the phenomenological model of pure silica, we may argue that the quantity da/drj may change 
with Na content because of different possible effects: i) a change of the elastic response of the system (via dO/drj), ii) 
a change in the functional dependence of the BO polarizability on the SiOSi angle (cfr. Fig. 8), iii) other structural 
parameters which become relevant with the presence of Na (such as the NaO interaction for instance) control the 
response to strain of the polarizability of BO, iv) the NBO contribute to a term in da/drj. In ref.^ Lines has con- 
jectured that the decrease in the photoelastic coefficients in NS glasses might be due to an increased modulation of 
the SiO bond length by strain, i.e. the Na ions would simply modify the mechanical response of the glass. However, 
we can recognize that this conjecture is not supported by the results in table I. In fact, the modulation of the SiO 
bond length upon even decreases in NS3 with respect to pure silica (cfr. table I). In order to clarify these issues, we 
have extended the phenomenological model for pure silica by considering different polarizabilitics for NBO and BO. 
We have introduced an additional parameter, aNBO, which describes a spherical polarizability of NBO ions. We have 
further assumed that aNBO may depend on the local electric field (Eioc) produced by the Na+ ions. The modulation 
of the Na-NBO distances upon strain would thus contribute to da/drj in the photoelastic tensor via a term of the 
form 

duNBO _ duNBO dEioc , , 

dr^ ~ dEioe drj ' ^ ' 

We have fitted the parameter daNBo/dEioc, which describes the response of the NBO polarizability to the local field 
of Na+ on the photoelastic coefficients of a sodium silicate crystal: the natrosilite Na2Si205^^. 



1. Natrosilite 



Natrosilite is a layered material with NBO nearly aligned along the c-axis, perpendicular to the siloxane layers (Fig 
1). Details on the structure of natrosilite are given in the appendix. 




FIG. 9. Side view of the natrosilite crystal. Large (small) dark grey spheres are Na (O) ions. Light grey spheres are Si ions. 

By changing the length of the c axis with the other lattice parameters (a, b and /3, see appendix) fixed, the SiOSi 
angles of BO do not change while there is a large change of the Na-NBO distances. We have thus assumed that the 
change in the dielectric constant of natrosilite upon the ij^ strain would be entirely due to the change of aj^BO which 
would obviuosly lead to an overestimation of daNBo/dr]. We have thus considered a model for the dielectric response 
of natrosilite in which the BO have the same polarizability as in pure silica (cfr. section Va) and the two additional 
parameters of NBO, a^BO and da^Bo/dris, have been fitted on the ab-initio photoelastic coefficient of natrosilite 
Pi3: P23, and P33. We obtain a^BO = 14.6 a.u. and daNBo/dV3 = 20.5 a.u. From Eq. 11 
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where we have assumed 
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Npfa is the number of Na ions nearest neighbor to NBO and rpfaO is the Na-NBO distance. The average is over 
the nearest neighbor ions. In natrosihte NNa='i and < r^y^Q >= 0.177 which finaUy yields daNBo/dEioc = 
—29.3 a.u.. Table IV reports the ab-initio photoelastic coefRcients of natrosilite and the results of the phenomenological 
model outlined above. The quantitity daNBo/dEioc is negative which means that by moving the Na ions further 
away from the NBO its polarizability increases. In fact, by decreasing the local electric field on NBO its charge 
would become more diffuse and thus more polarizable. The ab-initio photoelastic tensor has been computed by finite 
differences from the dielectric tensor calculated within DFPT with the codes PWSCF and PHONONS as described 
in section II. The calculations have been performed at the experimental equilibrium volume (see appendix) . 

TABLE IV. Average values of the phenomenological polarizability of BO and NBO ions and photoelastic coefficients of 
natrosilite computed ab-initio (DFPT) and with the phenomenologial model described in the text. 



Model DFPT 



obo 10.6 

Wnbo 14.62 

£11 2.454 2.412 

£22 2.467 2.443 

£33 2.28 2 2.361 

pia 0.140 0.131 

P23 0.143 0.138 

P33 0.08 7 0.104 
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2. NS3 glass 



In the development of a model for photoclasticity in NS3, we have used the value of da^ bo / dEioc obtamed from 
the fitting on natrosihte described in the previous section. The differences in the environment of NBO and in the 
structural response to strain of NS3 with respect to natrosihte is accounted for by the term dEioc/dr] in Eq. 12. In 
our model of NS3, N^a—'^-S' and < f^^^ >= 0.166 (cfr. section III). To extend the phenomenological model to 
NS3, we have considered two extreme cases as described below. 

a) In the first case (model A) we have assumed that the polarizability of BO ions {aso) is the same as in pure 
a-Si02 and assigned by the curves in Fig. 8. The polarizability of NBO, ajyfBO, is then assigned by fitting the ab-initio 
dielectric constant of NS3 within our extended phenomenologial model. The contribution from daNBo/dEioc, fitted 
on natrosihte, is added. The results are reported in table V. The fitting yields a^BO — 17.0 a. u., a value larger than 
otNBO in natrosihte (cfr. table IV), as we would have expected by considering that daNBO / dEioc is negative and 
Eioc is lower in NS3 with respect to natrosihte (cfr. Eq. 13). If we change by 20 % the value of a^BO (a change 
comparable with the difference in an bo between natrosihte and NS3) in the fitting procedure on natrosihte, one 
obtains a value for daNBo/dEioc with a similar change of 20 % with respect to the data in table IV. However, a 
change in da^ bo /dEioc of the order of 20 % does not affect the the results on the photoelastic coefficients of NS3 
within the figures reported in table V. 

b) In the second case (model B) ajsiBO is set to zero and aBO is still assigned by the function given in Fig. 8, 
but rescaled by a multiplicative factor in order to reproduce the dielectric constant of NS3. In this way we have also 
rescaled the term daBo/drj depending on the SiOSi angle. The increase of the polarizability of the BO upon Na 
insertion is supported by the comparison of the calculated Born effective charges for our models of NS3 and a-Si02 
shown in Fig. 10. The presence of Na induces a nearly uniform increase of the Born effective charges of BO ions. The 
valence electrons of the ionized Na atoms are thus transferred to both NBO and BO ions. A larger charge on the BO 
ions implies a larger polarizability. 

-1.3 
-1.4 
-1.5 

A 

1^-1.6 

V 

-1.7 
-1.8 



110 120 130 140 150 160 170 180 
Si-O-Si angle [degree] 

FIG. 10. Dependence of the Born effective charges of oxygen atoms as a function of the Si-O-Si angle in the models of NS3 
(black circles) and a-Si02 (red triangles, from ref.®). BO and NBO indicate bridging and non-bridging oxygen ions, respectively. 

In this model the contribution from NBO is added as well with the same parameter da n bo /dEioc fitted on natrosihte 
and used in model A. The results are reported in table V. 

TABLE V. Average values of the phenomenological polarizability of BO and NBO ions and photoelastic coefficients according 
to the two models (A and B) described in the text. The ab-initio (DFPT) photoelastic coefficients of NS3 and pure a-Si02 
(ref.^) are also reported. Results including and neglecting the term doNBO /dEioc (in a.u., see text) are both reported. 








Model A 


NS3 


Model B 


DFPT 


a-Si02 
Model 


DFPT 


Qbo 


10.6 


10.6 


17.8 


17.8 




10.6 




OINBO 


17.0 


17.0 














doLN BO /dEioc 





-29.33 





-29.33 








P21 


0.25 


0.24 


0.24 


0.23 


0.167 


0.227 


0.220 


P22 


0.18 


0.17 


0.14 


0.13 


0.072 


0.097 


0.057 
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Both models A and B largely overestimate the photoelastic coefficients in NS3 with respect to pure a-Si02. 

As mentioned before the value used for daNso/dEioc overestimates the real contribution of NBO to da/dt] in NS3. 
Yet, although overestimated, NBO gives a very small contribution to the photoelastic tensor. Even within model B 
where the contribution of BO is also rescaled by a factor 1.7 by still keeping aso dependent on the SiOSi angle only, 
the model photoelastic coefficients are sizably larger than the ab-initio data. The failure of these models suggests that 
either the shape of the functions in Fig. 8 changes upon Na insertion and/or other structural parameters, in addition 
to the SiOSi angle, influence the polarizability of BO ions. 

VI. CONCLUSIONS 

Photoelasticity in a model of sodium silicate glass (NS3) has been studied within density functional perturbation 
theory. The NS3 model is generated by quenching from the melt in combined classical and Car-Parrinello molecular 
dynamics simulations. 

The calculated photoelastic coefficients are in good agreement with experimental data and further confirm the 
reliability of DFPT already assessed for pure silica polymorphs in our previous work®. In particular, the calculation 
reproduces quantitatively the decrease of the photoelastic response induced by the insertion of Na, as measured 
experimentally''. Aiming at identifying the microscopic mechanisms thorough which sodium modifies the photoelastic 
response of the glass, we have extended to NS3 a phenomenological model of photoelasticity developed for pure silica 
in our previous work®. It comes out that the contribution of NBO to the photoelastic tensor (via a modulation of the 
NBO polarizability with strain) is not large enough to explain the decrease of the photoelastic coefficient of NS3 with 
respect to that of pure a-Si02. Moreover, although a charge transfer takes place from the ionized Na ions to the BO 
ions, a simple increase of the polarizability of BO is not sufficient to explain the change in the photoelastic response 
by keeping the SiOSi angle as the only structural parameter which modules the polarizability upon strain as in pure 
a-Si02. The modulation upon strain of other structural parameters should be called for. In order to quantify these 
latter effects, a more detailed modeling of the must be devised, requiring the fitting of the dielectric properties over 
a large database of sodium silicate crystals. 

* present address: Computational Science, Department of Chemistry and Applied Biosciences, ETH Zurich, USI 
Campus, Via Giuseppe Buffi 13, CH-6900 Lugano, Switzerland. 
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APPENDIX A 

Natrosilite, Na20(Si02)2, is a phyllosilicate crystal with space group P21/c and 4 formula units per unit celP®. 
We have optimized the internal structure at the experimental equilibrium lattice parameters a=8.13 A, 6=4.85 A, 
c=12.33 A and /? = 104.3°26. We have used the code PWSCF^* and 2x2x2 Monkhorst-Pack^^ mesh for Brillouin 
Zone integration. The experimental and theoretical positions (in crystallographic units) of the 9 independent atoms 
are reported in table VI 



TABLE VI. Experimental^® and theoretical positions (in crystallographic units) of the independent atoms of natrosilite. 





X 


Exp. 

y 


z 


X 


DFT 

y 


z 


Si 


0.028 


0.184 


0.183 


0.027 


0.164 


0.181 


Si 


0.403 


0.295 


0.277 


0.402 


0.304 


0.276 


Na 


0.379 


0.753 


0.443 


0.381 


0.752 


0.443 


Na 


0.137 


0.225 


0.473 


0.143 


0.210 


0.475 


O 


0.029 


0.859 


0.215 


0.026 


0.832 


0.216 


o 


0.454 


0.620 


0.267 


0.445 


0.637 


0.262 





0.226 


0.246 


0.181 


0.225 


0.242 


0.179 


o 


0.391 


0.232 


0.401 


0.388 


0.243 


0.400 


o 


0.093 


0.755 


0.436 


0.093 


0.721 


0.439 
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